
import geopandas as pd
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import imageio # Used for making animated GIFs
import numpy as np
from IPython.display import Image
from osgeo import gdal # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.warp
plt.rcParams['figure.figsize'] = (20, 20)
os.listdir("input")

/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(

['lds-nz-road-centrelines-topo-150k-FGDB.zip',
 'lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip',
 'statsnzregional-council-2021-clipped-generalised-FGDB.zip',
 'statsnz2018-census-usually-resident-population-and-age-groups-by-te-FGDB.zip',
 'Subnational population projections 2018-2048.xls',
 'lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip']


%%time
REGC = pd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)

epsg:2193
CPU times: user 1.75 s, sys: 986 ms, total: 2.73 s
Wall time: 14.6 s


%%time
df = pd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")

CPU times: user 1min 46s, sys: 3.2 s, total: 1min 49s
Wall time: 1min 49s


print(df.columns)
print(df.crs)
display(df.sample(5))

Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
       'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
       'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
       'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
       'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
      dtype='object')
epsg:2193


%%time
df = pd.clip(df, AKL)

CPU times: user 46.4 s, sys: 8.55 ms, total: 46.4 s
Wall time: 46.4 s


df.sample(5)


df.Name_2018.value_counts()

Exotic Forest                                3981
Indigenous Forest                            3673
Manuka and/or Kanuka                         2282
Broadleaved Indigenous Hardwoods             1788
Built-up Area (settlement)                   1350
High Producing Exotic Grassland              1326
Mangrove                                     1151
Urban Parkland/Open Space                    1099
Estuarine Open Water                          441
Orchard, Vineyard or Other Perennial Crop     436
Short-rotation Cropland                       362
Lake or Pond                                  326
Herbaceous Saline Vegetation                  303
Low Producing Grassland                       291
Gorse and/or Broom                            287
Forest - Harvested                            266
Sand or Gravel                                252
Deciduous Hardwoods                           201
Surface Mine or Dump                          132
Mixed Exotic Shrubland                        120
Herbaceous Freshwater Vegetation              118
Transport Infrastructure                      107
River                                          15
Flaxland                                        9
Gravel or Rock                                  9
Fernland                                        2
Matagouri or Grey Scrub                         1
Name: Name_2018, dtype: int64


def simplify_classes(code):
    if code in [1, 2, 5]:
        return 1, "Urban"
    elif code in [68,69,71]:
        return 2, "Vegetation"
    elif code in [0,20,21,22,45,46]:
        return 3, "Water"
    else:
        return 4, "Other"

summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
    print(year)
    class_year = f"Class_{year}"
    df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
    df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
    summary.append(df[class_year + "_simplified_name"].value_counts())

1996
2001
2008
2012
2018


pd.GeoDataFrame(summary).plot.area()

<AxesSubplot:>


%%capture
# %%capture suppresses output
if not os.path.isfile("land_use.gif"):
    ims = []
    years = [1996, 2001, 2008, 2012, 2018]
    for year in years:
        ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
        ax.set_title(year)
        ax.figure.tight_layout()
        canvas = ax.figure.canvas
        canvas.draw() # draw the canvas, cache the renderer
        image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
        image = image.reshape(canvas.get_width_height()[::-1] + (3,))
        ims.append(image)
    imageio.mimsave("land_use.gif", ims, fps=1)


with open('land_use.gif','rb') as file:
    display(Image(file.read()))


cols = [f"Class_{year}_simplified_code" for year in years]
cols

['Class_1996_simplified_code',
 'Class_2001_simplified_code',
 'Class_2008_simplified_code',
 'Class_2012_simplified_code',
 'Class_2018_simplified_code']


%%time
geocube = make_geocube(
    vector_data=df,
    output_crs="epsg:2193",
    measurements=cols,
    resolution=(-100, 100),
    fill=0, # NaNs, like offshore areas, will be 0
)
geocube

CPU times: user 18.4 s, sys: 9.83 ms, total: 18.4 s
Wall time: 18.5 s

<xarray.Dataset>
Dimensions:                     (x: 1001, y: 1320)
Coordinates:
  * y                           (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                           (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                 int64 0
Data variables:
    Class_1996_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2001_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2008_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2012_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2018_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


geocube.Class_2018_simplified_code.plot()

<matplotlib.collections.QuadMesh at 0x7f6c907ea910>


for year in years:
    print(year)
    outfile = f"output/land_use_{year}.tif"
    if not os.path.isfile(outfile):
        geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize

1996
2001
2008
2012
2018


%%time
pop = pd.read_file("input/statsnz2018-census-usually-resident-population-and-age-groups-by-te-FGDB.zip!2018-census-usually-resident-population-and-age-groups-by-te.gdb")
pop = pd.clip(pop, AKL)
pop.head()

CPU times: user 4.16 s, sys: 0 ns, total: 4.16 s
Wall time: 4.15 s


pop.columns

Index(['TALB2018_V1_00', 'TALB2018_V1_00_NAME', 'LAND_AREA_SQ_KM',
       'AREA_SQ_KM', 'Age_life_cycle_groups_15_29_years_2006',
       'Age_life_cycle_groups_30_64_years_2006',
       'Age_life_cycle_groups_65_years_and_over_2006',
       'Age_life_cycle_groups_total_2006',
       'Age_life_cycle_groups_Under_15_years_2006',
       'Age_life_cycle_groups_15_29_years_2013',
       'Age_life_cycle_groups_30_64_years_2013',
       'Age_life_cycle_groups_65_years_and_over_2013',
       'Age_life_cycle_groups_total_2013',
       'Age_life_cycle_groups_Under_15_years_2013',
       'Age_life_cycle_groups_15_29_years_2018',
       'Age_life_cycle_groups_30_64_years_2018',
       'Age_life_cycle_groups_65_years_and_over_2018',
       'Age_life_cycle_groups_total_2018',
       'Age_life_cycle_groups_Under_15_years_2018', 'Pop_Total_2006',
       'Pop_Total_2013', 'Pop_Total_2018',
       'Age_life_cycle_groups_15_29_years_2013_2018_percentChange',
       'Age_life_cycle_groups_30_64_years_2013_2018_percentChange',
       'Age_life_cycle_groups_65_years_and_over_2013_2018_percentChan',
       'Age_life_cycle_groups_total_2013_2018_percentChange',
       'Age_life_cycle_groups_Under_15_years_2013_2018_percentChange',
       'Shape_Length', 'geometry'],
      dtype='object')


cols = ["Pop_Total_2006", "Pop_Total_2013", "Pop_Total_2018"]
d = pop[cols].to_numpy()
d.min(), d.mean(), d.max()

(39, 63916.25, 140970)


%%capture
# %%capture suppresses output
if not os.path.isfile("pop.gif"):
    ims = []
    for col in cols:
        ax = pop.plot(column=col, legend=True, vmin=0, vmax=d.max(), cmap='OrRd')
        ax.set_title(col)
        ax.figure.tight_layout()
        canvas = ax.figure.canvas
        canvas.draw() # draw the canvas, cache the renderer
        image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
        image = image.reshape(canvas.get_width_height()[::-1] + (3,))
        ims.append(image)
    imageio.mimsave("pop.gif", ims, fps=1)


with open('pop.gif','rb') as file:
    display(Image(file.read()))


%%time
popcube = make_geocube(
    vector_data=pop,
    measurements=["Pop_Total_2006", "Pop_Total_2013", "Pop_Total_2018"],
    like=geocube, # Ensures dimensions match
    fill=0 # NaNs, like offshore areas, will be 0
)
popcube

CPU times: user 271 ms, sys: 19.5 ms, total: 290 ms
Wall time: 289 ms

<xarray.Dataset>
Dimensions:         (x: 1001, y: 1320)
Coordinates:
  * y               (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x               (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref     int64 0
Data variables:
    Pop_Total_2006  (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    Pop_Total_2013  (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    Pop_Total_2018  (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


for col in cols:
    outfile = f"output/{col}.tif"
    if not os.path.isfile(outfile):
        # byte max value is 255, and we have larger values than that here. uint32 max value is 4294967295
        popcube[col].rio.to_raster(outfile, dtype=np.uint32)


%%time
roads = pd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")

CPU times: user 8.6 s, sys: 20 ms, total: 8.62 s
Wall time: 8.62 s


%%time
akl_roads = pd.clip(roads, AKL)

CPU times: user 25.2 s, sys: 0 ns, total: 25.2 s
Wall time: 25.2 s


# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway


mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()


mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()


ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)


%%time
mway_cube = make_geocube(
    vector_data=mway,
    measurements=["lane_count"],
    like=geocube, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube

CPU times: user 259 ms, sys: 0 ns, total: 259 ms
Wall time: 256 ms

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    lane_count   (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
    mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)


src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


%%time
cbd = pop2018[pop2018.MB2020_V2_00 == "0433501"].copy()
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
    vector_data=cbd,
    like=geocube, # Ensures dimensions match
    fill=0
)
outfile = "output/cbd.tif"
if not os.path.isfile(outfile):
    cbd_cube.General_Electoral_Population.rio.to_raster(outfile, dtype=np.byte)
    
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)
CPU times: user 441 ms, sys: 10 ms, total: 451 ms
Wall time: 446 ms

Text(0.5, 1.0, 'Distance (m)')


bounds = AKL.total_bounds.tolist()
bounds

[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]


zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles

['EJ.tif',
 'DM.tif',
 'EL.tif',
 'DL.tif',
 'DJ.tif',
 'FK.tif',
 'DK.tif',
 'EK.tif',
 'FL.tif']


tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
DEM, transform = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100), dtype=np.int16)


print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(np.where(DEM>=0, DEM, np.nan), cmap='terrain')

-32767 -18316.517014943143 697 (1, 1320, 1018)

<AxesSubplot:>


cbd_dist = rasterio.open("output/cbd_dist.tif")
cbd_dist.read().shape

(1, 1320, 1001)


NODATA = -32767
warped_DEM = np.full_like(cbd_dist.read(), NODATA, dtype=np.int16)

rasterio.warp.reproject(DEM,
                        warped_DEM,
                        src_transform=transform,
                        src_nodata=NODATA,
                        dst_nodata=NODATA,
                        src_crs=tile_datasets[0].crs,
                        dst_crs=cbd_dist.crs,
                        dst_transform=cbd_dist.transform
)
print(DEM.shape, warped_DEM.shape)

(1, 1320, 1018) (1, 1320, 1001)


meta = tile_datasets[0].meta
print(meta)
meta.update({
    "dtype": "int16",
    "height": cbd_dist.height,
    "width": cbd_dist.width,
    "transform": cbd_dist.transform
})
print(meta)
outfile = "output/dem.tif"
if not os.path.isfile(outfile):
    with rasterio.open(outfile, "w", **meta) as dest:
       dest.write(warped_DEM)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
       0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703900.0,
       0.0, -100.0, 6002400.0)}


gdal.DEMProcessing('output/slope.tif', outfile, 'slope')
slope = rasterio.open("output/slope.tif").read()
print(np.nanmin(slope), np.nanmean(slope), np.nanmax(slope), slope.shape)
rasterio.plot.show(np.where(slope>=0, slope, np.nan), cmap='terrain')

-9999.0 -5778.3354 52.492313 (1, 1320, 1001)

<AxesSubplot:>


rasters = [rasterio.open(f"output/{f}") for f in os.listdir("output") if f.endswith(".tif")]
params = set([(r.shape, r.res, r.crs, r.count, r.bounds) for r in rasters])
print(params)
# Assert all rasters have the same shape, pixel size, CRS, number of bands, and bounds
assert len(params) == 1

{((1320, 1001), (100.0, 100.0), CRS.from_epsg(2193), 1, BoundingBox(left=1703900.0, bottom=5870400.0, right=1804000.0, top=6002400.0))}


!ls -Ggh output

total 27M
-rw-r--r-- 1 1.3M Apr 16 14:44 cbd.tif
-rw-r--r-- 1 2.6M Apr 16 16:09 cbd_dist.tif
-rw-r--r-- 1 2.6M Apr 16 16:46 dem.tif
-rw-r--r-- 1 1.3M Apr 16 09:50 land_use_1996.tif
-rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2001.tif
-rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2008.tif
-rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2012.tif
-rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2018.tif
-rw-r--r-- 1 1.3M Apr 16 09:50 mway.tif
-rw-r--r-- 1 2.6M Apr 16 16:09 mway_dist.tif
-rw-r--r-- 1 2.6M Apr 16 09:50 pop2013.tif
-rw-r--r-- 1 2.6M Apr 16 09:50 pop2018.tif
-rw-r--r-- 1 5.1M Apr 16 16:52 slope.tif

	TALB2018_V1_00	TALB2018_V1_00_NAME	LAND_AREA_SQ_KM	AREA_SQ_KM	Age_life_cycle_groups_15_29_years_2006	Age_life_cycle_groups_30_64_years_2006	Age_life_cycle_groups_65_years_and_over_2006	Age_life_cycle_groups_total_2006	Age_life_cycle_groups_Under_15_years_2006	Age_life_cycle_groups_15_29_years_2013	...	Pop_Total_2006	Pop_Total_2013	Pop_Total_2018	Age_life_cycle_groups_15_29_years_2013_2018_percentChange	Age_life_cycle_groups_30_64_years_2013_2018_percentChange	Age_life_cycle_groups_65_years_and_over_2013_2018_percentChan	Age_life_cycle_groups_total_2013_2018_percentChange	Age_life_cycle_groups_Under_15_years_2013_2018_percentChange	Shape_Length	geometry
2	00300	Kaipara District	3108.710953	3114.579850	2631	8568	2715	18135	4215	2670	...	18135	18963	22869	25.6	15.2	33.8	20.6	16.4	902974.494161	MULTIPOLYGON (((1718576.593 5980699.195, 17185...
4	01200	Hauraki District	1270.087068	1270.087068	2568	8274	3051	17856	3963	2634	...	17856	17808	20022	12.3	10.2	21.3	12.4	7.5	311223.430113	POLYGON ((1801726.041 5899678.947, 1803788.474...
5	01300	Waikato District	4403.282847	4516.062518	9888	27450	5580	57588	14670	10695	...	57588	63378	75618	23.8	17.2	26.9	19.3	16.6	491881.373163	POLYGON ((1747419.052 5870470.171, 1747413.783...
29	07606	Henderson-Massey Local Board Area	53.219988	64.608373	21441	43914	9264	98790	24174	23274	...	98790	107685	118422	12.7	10.6	8.3	10.0	7.0	47587.459227	MULTIPOLYGON (((1748774.713 5920936.996, 17483...
30	07607	Waitakere Ranges Local Board Area	305.732502	943.004892	8343	22752	3372	45498	11031	8763	...	45498	48399	52095	13.4	5.6	18.6	7.6	3.0	148724.962888	POLYGON ((1742888.086 5900837.328, 1742806.519...

	t50_fid	name_ascii	macronated	name	hway_num	rna_sufi	lane_count	way_count	status	surface	geometry
512	100120610	KAIPARA COAST HIGHWAY	N	KAIPARA COAST HIGHWAY	16	3007739	2	None	None	sealed	LINESTRING (1732000.000 5944172.070, 1732048.5...
2933	3198057	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748581.508 5968975.145, 1748558.4...
2934	3198059	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748171.047 5971284.152, 1748129.9...
3320	3200754	PAERATA ROAD	N	PAERATA ROAD	22	3000260	2	None	None	sealed	LINESTRING (1767236.112 5888088.508, 1767244.3...
3324	3200792	UPPER HARBOUR MOTORWAY	N	UPPER HARBOUR MOTORWAY	18	3047073	4	None	None	sealed	LINESTRING (1747954.314 5927269.837, 1747970.0...
...	...	...	...	...	...	...	...	...	...	...	...
138240	100048291	AUCKLAND-WAIWERA MOTORWAY	N	AUCKLAND-WAIWERA MOTORWAY	1	3067966	7	None	None	sealed	LINESTRING (1755881.018 5922863.734, 1755886.4...
138301	100048432	AUCKLAND-HAMILTON MOTORWAY	N	AUCKLAND-HAMILTON MOTORWAY	1	3017109	1	None	None	sealed	LINESTRING (1765115.647 5909916.697, 1765092.7...
138337	100048532	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	4	None	None	sealed	LINESTRING (1748892.089 5949596.727, 1748892.0...
138369	100048589	PORT ALBERT ROAD	N	PORT ALBERT ROAD	16	3013274	2	None	None	sealed	LINESTRING (1734173.019 5980575.187, 1734175.6...
138680	100118365	SOUTH-WESTERN MOTORWAY	N	SOUTH-WESTERN MOTORWAY	20	3018532	4	None	None	sealed	LINESTRING (1760066.252 5908184.133, 1760043.9...

Land use¶

Population density¶

Distance from major roads¶

Distance from CBD¶

Slope¶

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
174810	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	69	69	69	69	69	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000180191	MULTIPOLYGON (((1271800.779 4844357.101, 12717...
4423	Gorse and/or Broom	Gorse and/or Broom	Gorse and/or Broom	Gorse and/or Broom	Gorse and/or Broom	51	51	51	51	51	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000105107	MULTIPOLYGON (((1528090.924 5194424.609, 15280...
430719	Lake or Pond	Lake or Pond	Lake or Pond	Lake or Pond	Lake or Pond	20	20	20	20	20	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000045836	MULTIPOLYGON (((1261101.499 4966013.291, 12611...
461252	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	69	69	69	69	69	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000169706	MULTIPOLYGON (((1710651.930 6020180.126, 17106...
344602	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	yes	yes	yes	yes	yes	yes	Landcare Research	2014-06-30T00:00:00	lcdb1000098963	MULTIPOLYGON (((1612278.405 6152038.608, 16122...

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
173780	Urban Parkland/Open Space	Urban Parkland/Open Space	Urban Parkland/Open Space	Urban Parkland/Open Space	Urban Parkland/Open Space	2	2	2	2	2	...	no	yes	yes	yes	yes	yes	Landcare Research	2004-06-30T00:00:00	lcdb1000006473	POLYGON ((1757451.024 5944934.676, 1757463.730...
57529	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000093800	POLYGON ((1733473.489 5921727.296, 1733459.367...
408632	Low Producing Grassland	Low Producing Grassland	Exotic Forest	Exotic Forest	Exotic Forest	41	41	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2019-12-01T00:00:00	lcdb1000415656	POLYGON ((1778637.742 5904919.429, 1778637.742...
394053	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000094687	POLYGON ((1735422.343 5948607.284, 1735382.309...
507426	Herbaceous Saline Vegetation	Herbaceous Saline Vegetation	Herbaceous Saline Vegetation	Herbaceous Saline Vegetation	Herbaceous Saline Vegetation	46	46	46	46	46	...	yes	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000066435	POLYGON ((1749641.378 5917711.857, 1749687.537...